Métodos Estadísticos

# Población: ventas diarias ~ Uniforme(0, 20), media = 10
set.seed(9915)
tibble(x = runif(n=1e9, min=0, max=20)) |> pull(x) |> mean()

# Muestra de solo 25 días
set.seed(9915)
tibble(x = runif(n=25, min=0, max=20)) |> pull(x) |> mean()

Regresión Lineal: Estimación

• Objetivo: estimar los coeficientes poblacionales desconocidos a partir de una muestra

• Los coeficientes estimados minimizan la Suma Cuadrática de Residuos, \(\small \hat{e}\) (errores de estimación): es decir , las distancias entre

  • los datos observados: \(\small y_i\)

  • y los datos predichos por el modelo estimado \(\small \hat{y}_i=\hat{\beta}_0+\hat{\beta}_1 X_1+ \dots + \hat{\beta}_k X_p\)

• Esto equivale a las condiciones derivadas de suponer \(\small E(\varepsilon|X)=0\)

Regresión Lineal: Estimación (cont.)

• Los parámetros del modelo \(f(.)\) estimado, los \(\small \beta_j\), son estadísticos con variabilidad muestral, medida por su error estándar \(\small se(\widehat{\beta}_j)\)

  • Usados para construir intervalos de confianza y para contrastar hipótesis sobre los parámetros poblacionales, p.e., significatividad (\(\small \beta_1=0\)) y signo.

• La función lm() con la notación de fórmula crea un objeto lista con resultados

descuento <- read.csv("data/descuento.csv") |>
        mutate(zona = parse_factor(zona), mujer = as.factor(mujer))
modelo <- lm(data = descuento, ventas ~ renta + descuento)
(sum.modelo <- summary(modelo))
nobs(modelo)

• La predicción \(\widehat{y}\) también está sujeta a incertidumbre por la estimación: se puede calcular su error estándar e intervalos de confianzas

pred <- predict(modelo, type = "response",
                se.fit = T, interval = "confidence")
cbind(descuento$ventas, pred$fit, pred$se.fit) |> head()

Ampliando el Modelo de Regresión Lineal

• Incorporar información cualitativa: incluimos una variable binaria (dummy) para cada categoría, excepto una (grupo de referencia)

  • Cada coeficiente es el efecto diferencial de ese grupo respecto al grupo de referencia

  • Las dummies se crean automáticamente con factores (sin orden)

    • también se pueden usar la biblioteca fastDummies

lm(data = descuento, ventas ~ renta + descuento + zona) |> summary()

Regresión Lineal: comentarios finales

• Los resultados de (todos) estos modelos informan de la relación entre los regresores y la variables dependiente

• “Correlación no implica causalidad”, salvo en ensayos científicos aleatorios controlados, también llamados pruebas A/B

lm(data = descuento, ventas ~ descuento) |> summary()
lm(data = descuento, ventas ~ descuento + renta) |> summary()
lm(data = descuento, ventas ~ renta) |> summary()
lm(data = descuento, ventas ~ renta + educ) |> summary()

Regresión Logística

• Un modelo de regresión lineal para respuestas binarias no es adecuado porque genera predicciones fuera del rango \(\small [0,1]\)

• De manera que \(\small \Pr(Y=1|X)= f(x)= \color{blue}{\Lambda( \beta_0 + \beta_1 x_1 + \ldots + \beta_k x_k)} \in [0,1]\)

Regresión Logística (cont.)

• La regresión logística se puede estimar en R con la función glm() (similar a lm()); podemos incluir variables explicativas tanto cuantitativas como cualitativas, interacciones, etc.

censo <- read_csv("data/census.csv") |>
  mutate(across(where(is.character), ~parse_factor(.x)))

logit <- glm(data = censo, income ~  race + age + I(age^2) +
               log(hours_per_week)*sex +  education_1  +
               capital_gain + capital_loss, family = "binomial")
logit |>  summary()

• Como en el modelo lineal, podemos calcular predicciones de la probabilidad tanto en los mismos datos como en unos nuevos

predict(logit, type="response")

logit2 <- logit <- glm(data = censo, income ~  age + education_1,
         family = "binomial")
summary(logit2)
predict(logit2, type="response",
        newdata = tibble(age=c(20,60), education_1=c(9,16)))

Regresión Logística con más de dos clases

• La regresión logística se puede generalizar a situaciones con múltiples clases (modelos multinomiales) con un índice lineal para cada clase \[ \small \Pr(Y=c|X)=\frac{e^{\beta_{0c}+\beta_{1c}X_1+\dots+\beta_{kc}X_k}}{\sum_{l=1}^{C}e^{\beta_{0l}+\beta_{1l}X_1+\dots+\beta_{kl}X_k}} \]

• La librería glmnet() permite la estimación de estos modelos

health <- read_csv("data/health_insurance.csv") |>
  mutate(across(where(is.character), ~parse_factor(.x)))

library(glmnet)
x <- model.matrix(product ~ age + gender, data = health)
mod.glmnet <- glmnet(x = x, y = health$product, family = "multinomial",
                     lambda = 0, type.multinomial = "grouped")
coef(mod.glmnet)
predict(mod.glmnet, newx=x, type = "response")  # probabilidad de cada clase
predict(mod.glmnet, newx=x, type = "class")     # clase

Error de predicción

• Un modelo es mejor si sus predicciones se ajusten mejor a las observaciones

• El error de predicción es \(y - \widehat{y} = f(X) - \widehat{f}(X) + \varepsilon\)

  • \(f - \widehat{f}\) = error reducible (eligiendo modelo)

  • \(\varepsilon\) = error irreducible (variables no observadas)

• La función de pérdida (o coste) evalúa cómo valoramos las desviaciones

Métricas de error de predicción (cuantitativa)

• Mean Square Error (Error Cuadrático Medio): \(\small MSE(y,\widehat{y})={\frac{1}{n}\sum_{i=1}^{n}\left(y-\widehat{y}\right)^2}\)

  • penaliza grandes desviaciones

  • \(\small R^2\) y \(\small R^2\)-ajustado son variantes del MSE, pero solo sirven para comparar modelos con la misma variable dependiente.

• Root Mean Square Error: \(\small RMSE(y,\widehat{y})=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y-\widehat{y}\right)^2}\)

  • mismas unidades que \(\small y\)

• Mean Absolute Error: \(\small MAE(y,\widehat{y})=\frac{1}{n}\sum_{i=1}^{n}\left|y-\widehat{y}\right|\)

• Otras medidas basadas en distintas funciones de pérdida, la verosimilitud del modelo (\(\small AIC\), \(\small BIC\)), etc

Seleccionar el mejor modelo

• ¿Podemos predecir el número de visitantes en función de la temperatura?

library(mosaicData)
data(RailTrail)
RailTrail |> ggplot(aes(x = hightemp, y = volume)) +
  geom_point() +  geom_smooth()

• ¿Cuál es el mejor modelo?

  • \(\small volume = \beta_0 + \beta_1 hightemp + \varepsilon\)

  • \(\small volume = \beta_0 + \beta_1 hightemp + \beta_2 hightemp^2 + \varepsilon\)

  • …

  • \(\small volume = \beta_0 + \beta_1 hightemp + \dots + \beta_{22} hightemp^{22} + \varepsilon\)

RailTrail |> ggplot(aes(x = hightemp, y = volume)) +
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ poly(x,22) ) +
  coord_cartesian(ylim = c(100,750))

Muestras de entrenamiento y de prueba

• Las métricas de error (ej., \(\small MSE\)) suelen calcularse para predicciones de los mismos datos usados para ajustar/estimar el modelo (in-sample prediction)

  • Esta muestra se denomina muestra de entrenamiento (training sample)

• PERO queremos saber qué tal se predicen casos nuevos (out-sample prediction)

• Usar las métricas en muestras de entrenamiento implica problemas de “overfitting”: sobreajuste a las características de la muestra concreta

  • Un modelo menos flexible podría tener menor error de predicción con casos nuevos

• Debemos calcular las métricas de error con observaciones que el modelo NO ha usado antes: muestra de prueba (test sample)

“Overfitting”

MSE en la muestra de prueba

\[ \small E\left[\left(y-\widehat{f}(x)\right)^2\right] = E\left[\left(f(x)-\widehat{f}(x) + \varepsilon + E\left[\widehat{f}(x)\right]-E\left[\widehat{f}(x)\right] \right)^2\right] = \]

\[ \small =\underbrace{\left[E\left(\widehat{f}(x)\right)-f(x)\right]^2}_{(1)} + \underbrace{E\left(\left[\widehat{f}(x)-E\left(\widehat{f}(x)\right)\right]^2\right)}_{(2)}+Var(\varepsilon) \]

• \(\small (1)=\left[Sesgo\left(\widehat{f}(x)\right)\right]^2\): error por supuestos erróneos en \(f\)

  • ajuste insuficiente (underfit) al perder relaciones relevantes entre \(X\) e \(Y\)

• \(\small (2)=Var\left(\widehat{f}(x)\right)\): sensibilidad a fluctuaciones en el entrenamiento

  • si el algoritmo modela puro ruido en entrenamiento, ajustará bien allí, pero predecirá mal casos nuevos (overfit)

“Trade-off” Varianza–Sesgo

Medir el Error en la Clasificación

• Los modelos de clasificación NO predicen directamente la categoría, sino la probabilidad de que una observación pertenezca a cada categoría

• La clase predicha será aquella con mayor probabilidad. En el caso binario, implicar superar el umbral de 0.5 (se deben probar varios valores)

censo <- read_csv("data/census.csv") |>
  mutate(income = parse_factor(income))
logit <- glm(data = censo, income ~ capital_gain,
             family = "binomial")
prob.predict <- predict(logit, type = "response")

umbral <- 0.5
cat.predict  <- if_else(prob.predict > umbral, 1, 0)
cbind(censo$income, cat.predict, prob.predict) |> head(10)

• Como no tiene sentido diferencia de clases (variables categóricas), NO se pueden calcular medidas como el MSE y otros relacionados

Matriz de Confusión

• Tabular categorías observadas frente a las categorías predichas

table(cat.predict, censo$income)

• Existen varias medidas derivadas de la matriz de confusión

Métricas con la matriz de confusión

• Tasa de observaciones correctamente clasificadas (exactitud o accuracy)

\[ \scriptsize ACCUR=\frac{VP+VN}{VP+FP+VN+FN} \]

• No es informativo cuando algunas clases son infrecuentes (datos desequilibrados)

  • si hay poco fraude/enfermos (ej., 5%), predecir que nunca hay fraude implica \(\scriptsize ACCUR=95\%\), PERO NO detecta fraude/enfermedad

• El estadístico Kappa (\(\small \kappa\)) es una medida similar, pero que ajusta por lo se esperaría solo por azar (corrigiendo en parte el desequilibrio entre clases).

Métricas con la matriz de confusión (cont.)

• La tasa de verdaderos positivos o sensibilidad (recall) es el porcentaje de verdaderos positivos sobre el total de positivos observados \[ \scriptsize TVP=SENSIT=\frac{VP}{VP+FN} \]

  • ej., tasa de fraude/enfermos existentes que se detectan correctamente

• La tasa de verdaderos negativos o especificidad es el porcentaje de verdaderos negativos sobre el total de negativos observados \[ \scriptsize TVN=ESPECIF=\frac{VN}{VN+FP} \]

  • ej., tasa de “otras” opciones que se clasifican correctamente

  • Tasa de falsos positivos: \(\scriptsize TFP = 1 - TVN = 1 - ESPECIF\)

Métricas con la matriz de confusión (y 3)

• La exactitud equilibrada (Balanced Accuracy) es una media de la sensibilidad y de la especificidad

• La precisión o valor de predicción positivo es la cantidad de verdaderos positivos sobre el total de positivos predichos

\[ \scriptsize PREC=\frac{VP}{VP+FP} \]

• La familia de medidas \(\small F_{\beta}\) es una ratio de la importancia ponderada de la sensibilidad y de la precisión: \(\scriptsize F_{\beta}=\frac{(1+\beta)^2 \times SENSIT \times PREC}{\beta^2 \times SENSIT + PREC}\)

  • Para \(\scriptsize \beta<1\), se da menos importancia a la sensibilidad: los falsos positivos se consideran más costosos

  • Para \(\scriptsize \beta>1\), los falsos negativos son más costosos y para \(\scriptsize \beta=1\) son igualmente costosos

Curva ROC (“Receiver Operating Characteristic”)

AUC (“area under the curve”)

• La AUC es el área bajo la curva ROC: ofrece una medida agregada de rendimiento entre 0 (todas las clasificaciones incorrectas) y 1 (todas correctas)

• Resume la curva ROC y permite comparar curvas que se cruzan

Extensiones. Métricas adicionales

• Con más de dos clases, se realiza un análisis AUC-ROC para cada categoría (frente a las demás) y se promedian (ej., ponderando por casos en cada clase)

• Con clases desequilibradas, se puede preferir en lugar de la ROC un gráfico de precisión frente sensibilidad (precision-recall) y su correspondiente AUC (PR-AUC)

• Existen múltiples funciones de pérdida (o coste de clasificación) posibles.

  • Las relacionadas con la curva de ganancia consideran el coste de alcanzar un cierto nivel de sensibilidad

  • Otras se basan en la función de verosimilud o la entropía como medidas de pérdida (ej. mean log loss)

Evaluación de Modelos: entrenamiento y prueba

• Para minimizar problemas de underfit y, sobre todo, de overfit, DEBEMOS dividir aleatoriamente el conjunto de datos en dos partes:

• ¿Por qué renunciar a parte de los datos si sabemos que un tamaño muestral grande es importante? Evaluar correctamente un modelo lo es mucho más

• La estimación del error en prueba puede ser volátil dependiendo de las observaciones incluidas en cada grupo

Evaluación de Modelos: Validación cruzada

• Para evitar que los datos sean sensibles a una partición concreta, se usa validación cruzada (cross-validation o rotation estimation)

• Se repite varias veces y de forma ordenada el proceso de remuestreo para la partición en grupos de entrenamiento y prueba (similar a bootstrap)

• Permite utilizar todas las observaciones de la muestra, tanto para estimar como para evaluar el modelo (aunque no a la vez)

Validación cruzada de K bloques

Validación cruzada aleatoria (RCV) y LOOCV